      subroutine nonresmat(o,Mbj2,d1,d2,d3,d4,d5,Cp,Cap,Rp,Rap,PE)
      implicit none
c DECLARATIONS of VARIABLES
c declaration of dimensions
      integer d1,d2,d3,d4,d5
Cf2py intent(in) d1,d2,d3,d4,d5
c declaration of boundaries for the matrix
      integer c1i,c1e,c2i,c2e,c3i,c3e
c
c declaration of the matrix
c  output matrix
      real*8 o(1:d1,1:d1)
Cf2py intent(out) o
Cf2py depend(d1) o
c  coefficient matrix (element wise squared)
      real*8 Mbj2(1:d2,1:d2)
Cf2py intent(in) Mbj2
Cf2py depend(d2) Mbj2

c Matrix related to Ci\Cj
      real*8 Ci(1:d2,1:d2),Cj(1:d2,1:d2)
Cf2py depend(d2) Ci,Cj
      real*8 Cp,Cap,Calpha,Cbeta
Cf2py intent(in) Cp,Cap
      real*8 Cm11(1:d3*d4,1:d3*d4),Cm21(1:d3*d4,1:d3*d4)
      real*8 Cm12(1:d2,1:d2),Cm22(1:d2,1:d2)
      real*8 Cv1(d3*d4),Cv2(d5),Cv3(d2)
Cf2py depend(d3,d4) Cm11,Cm21,Cv1
Cf2py depend(d2) Cm12,Cm22,Cv3
Cf2py depend(d5) Cv2

c Matrix related to Rj\Rk
      real*8 Rj(1:d2,1:d2),Rk(1:d5, 1:d2)
Cf2py depend(d2) Rj
Cf2py depend(d2,d5) Rk
      real*8 Rv1(d3*d4),Rv2(d5)
      real*8 Rv3(d2)
Cf2py depend(d3,d4) Rv1
Cf2py depend(d5) Rv2
Cf2py depend(d2) Rv3
      real*8 Rp,Rap
Cf2py intent(in) Rp,Rap

c Matrix related to PEi\PEk
      real*8 PEi(1:d2,1:d5),PEk(1:d5,1:d5)
Cf2py depend(d2,d5) PEi
Cf2py depend(d5) PEk
      real*8 PE
Cf2py intent(in) PE
      real*8 PEv1(1:d3*d4)
c      real*8 PEm1(1:d2,1:d5)

c Operator matrixw
      real*8 Sz(1:2,1:2),Sd(1:2,1:2),Id(1:d5,1:d5)

      integer i,i2,io,j,k,l,m
      real*8 ms

c DEFINITIONS and INITIALIZATIONS
c constants
c  dimensions
      c1i=1
      c1e=d2
      c2i=c1e+1
      c2e=c1e+d2
      c3i=c2e+1
      c3e=c2e+d5

c matrix
      o=0
      Ci=0
      Cj=0
      Rj=0
      Rk=0
      PEi=0
      PEk=0
      Sz=0
      Sd=0
      Id=0

cccccccccccccccccccccccccccccccccccccccccccc
c build_sz(Sz,2)
cccccccccccccccccccccccccccccccccccccccccccc
      Sz=0
      ms=(2-1)/2.0
      do k=1,2,1
      Sz(k,k)=ms
      ms=ms-1
      end do


cccccccccccccccccccccccccccccccccccccccccccc
c build_sd(Sd,2)
cccccccccccccccccccccccccccccccccccccccccccc
      Sd=0
      do k=1,2,1
      Sd(k,k)=1
      end do

cccccccccccccccccccccccccccccccccccccccccccc
c build_sd(Id,d5)
cccccccccccccccccccccccccccccccccccccccccccc
      Id=0
      do k=1,d5,1
      Id(k,k)=1
      end do

cccccccccccccccccccccccccccccccccccccccccccc
c
c      call build_ci(Ci,d2,d3,d4,d5,Cp,Cap)
c
cccccccccccccccccccccccccccccccccccccccccccc

      Calpha = (Cp - Cap) / 2.0
      Cbeta = (Cp + Cap) / 2.0
c matrix
      Cm11=0
      Cm12=0
      Cm21=0
      Cm22=0

cccccccccccccccccccccccccccccccccccccccccccc
c call k22(Cm11,Sz,2,2,Sz,2,2)
cccccccccccccccccccccccccccccccccccccccccccc
      Cm11(1,1)=1.0/4
      Cm11(2,2)=-1.0/4
      Cm11(3,3)=-1.0/4
      Cm11(4,4)=1.0/4

cccccccccccccccccccccccccccccccccccccccccccc
c call k22(Cm21,Sd,2,2,Sd,2,2)
cccccccccccccccccccccccccccccccccccccccccccc
      Cm21(1,1)=1
      Cm21(2,2)=1
      Cm21(3,3)=1
      Cm21(4,4)=1

cccccccccccccccccccccccccccccccccccccccccccc
c call k22(Cm12,Cm11,4,4,Id,d5,d5)
c call k22(Cm22,Cm21,4,4,Id,d5,d5)
cccccccccccccccccccccccccccccccccccccccccccc
      do i=1,4,1
       do j=1,d5,1
        io=j+(i-1)*d5
        Cm12(io,io)=Cm11(i,i)*Id(j,j)
        Cm22(io,io)=Cm21(i,i)*Id(j,j)
       end do
      end do

      Ci = Calpha *  Cm12 + Cbeta * Cm22


cccccccccccccccccccccccccccccccccccccccccccc
c
c      call build_cj(Cj,Mbj2,d2,d3,d4,d5,Cp,Cap)
c
cccccccccccccccccccccccccccccccccccccccccccc
      Cv2=1
      Cv3=0
      Cv1(1)=Cp
      Cv1(2)=Cap
      Cv1(3)=Cap
      Cv1(4)=Cp

cccccccccccccccccccccccccccccccccccccccccccc
c call k11(Cv3,Cv1,d3*d4,Cv2,d5)
cccccccccccccccccccccccccccccccccccccccccccc
      do i=1,d3*d4,1
       do j=1,d5,1
        Cv3(j+(i-1)*d5)=Cv1(i)*Cv2(j)
       end do
      end do

      do i=1,d2,1
       do j=1,d2,1
        Cj(i,j)=Cv3(j)*Mbj2(j,i)
       end do
      end do



cccccccccccccccccccccccccccccccccccccccccccc
c
c call build_rj(o,Mbj2,d2,d3,d4,d5,Rp,Rap)
c
cccccccccccccccccccccccccccccccccccccccccccc
      Rv1(1)=Rp
      Rv1(2)=Rap
      Rv1(3)=Rap
      Rv1(4)=Rp
      Rv3=0
cccccccccccccccccccccccccccccccccccccccccccc
c call k11(Rv3,Rv1,d3*d4,Rv2,d5)
cccccccccccccccccccccccccccccccccccccccccccc
      do i=1,d3*d4,1
       do j=1,d5,1
        Rv3(j+(i-1)*d5)=Rv1(i)
       end do
      end do

c      do i=1,d2,1
c       do j=1,d2,1
c        Rj(i,j)=Rv3(j)*Mbj2(j,i)
c       end do
c      end do
      do i=1,d2,1
       do j=1,d2,1
        Rj(i,i)=Rj(i,i)+Rv3(j)*Mbj2(j,i)
       end do
      end do




cccccccccccccccccccccccccccccccccccccccccccc
c
c call build_rk(o,Mbj2,d2,d3,d4,d5,Rp,Rap)
c
cccccccccccccccccccccccccccccccccccccccccccc
      Rv3=0

      do i=1,d5,1
       Rv2=0
       Rv2(i)=1
cccccccccccccccccccccccccccccccccccccccccccc
c call k11(Rv3,Rv1,d3*d4,Rv2,d5)
cccccccccccccccccccccccccccccccccccccccccccc
       do l=1,d3*d4,1
        do m=1,d5,1
         Rv3(m+(l-1)*d5)=Rv1(l)*Rv2(m)
        end do
       end do

       do j=1,d2,1
        do k=1,d2,1
         Rk(i,j)=Rk(i,j)+Rv3(k)*Mbj2(k,j)
        end do
       end do
      end do
c      print*,Rk

cccccccccccccccccccccccccccccccccccccccccccc
c
c call build_pei(o,d2,d3,d4,d5,PE)
c
cccccccccccccccccccccccccccccccccccccccccccc
      PEv1=1
cccccccccccccccccccccccccccccccccccccccccccc
c call k21c(PEm1,PEv1,d3*d4,Id,d5,d5)
cccccccccccccccccccccccccccccccccccccccccccc
      do i=1,d3*d4,1
       do i2=1,d5,1
c       j2=i2...
c        do j2=1,d5,1
         PEi(i2+(i-1)*d5,i2)=PEv1(i)*Id(i2,i2)*PE/d2
c        end do
       end do
      end do

cccccccccccccccccccccccccccccccccccccccccccc
c
c call build_pek(o,d5,PE)
c
cccccccccccccccccccccccccccccccccccccccccccc
      do i=1,d5,1
       PEk(i,i)=PE / d5
      end do
c
c
      o(c1i:c1e, c1i:c1e) = -Ci
      o(c2i:c2e, c1i:c1e) = Cj
      o(c2i:c2e, c2i:c2e) = -Rj
      o(c3i:c3e, c2i:c2e) = Rk
      o(c1i:c1e, c3i:c3e) = PEi
      o(c3i:c3e, c3i:c3e) = -PEk

      end subroutine nonresmat



      subroutine resmat(o_n,o_f,dtm,f,w,mwe,ev,d1,d2)
c
      implicit none
c DECLARATIONS of VARIABLES
c declaration of dimensions
      integer d1,d2
Cf2py intent(in) d1,d2
c
c declaration of the matrix
c  output matrix
      real*8 o_n(1:d1,1:d1),o_f(1:d1,1:d1)
Cf2py intent(out) o_n
Cf2py intent(in) o_f
Cf2py depend(d1) o_n,o_f
c  coefficient matrix (element wise squared)
      real*8 dtm(1:d2,1:d2)
Cf2py intent(in) dtm
Cf2py depend(d2) dtm
      real*8 ev(d2)
Cf2py intent(in) ev
Cf2py depend(d2) ev
c
c declaration of some constant
      real*8 f,w,mwe
Cf2py intent(in) f,w,mwe
      real*8 dts,gaussian,intensity,center
c dipolar_transition_strength
      integer i,j,im,jm
c
      o_n=0
      do i=1,d2,1
       do j=1,d2,1
        center = abs(ev(i) - ev(j))
        intensity = dtm(i,j)**2
        dts= gaussian(intensity,center,w,f)
        im=d2+i
        jm=d2+j

        o_n(im,im)=o_n(im,im)-mwe*dts
        o_n(jm,jm)=o_n(jm,jm)-mwe*dts
        o_n(im,jm)=o_n(im,jm)+mwe*dts
        o_n(jm,im)=o_n(jm,im)+mwe*dts
        dts=0
        intensity=0
       end do
      end do
      o_n=o_n+o_f
      end subroutine resmat



      function lorentzian(i,c,w,x)
      implicit none
c w the width
c c the center
      real*8 lorentzian
      real*8 i,c,w,x
Cf2py intent(in) i,c,w,x
      real*8 num,denom,pi
      num=0
      denom=0
      pi=3.1415926535
      num=i * (w / (2 * pi))
      denom=((x-c)**2+(w/2)**2)
      lorentzian=num/denom
      end function lorentzian
c
c
      function gaussian(i,c,w,x)
      implicit none
c w the width
c c the center
      real*8 gaussian,pi
      real*8 i,c,w,x
Cf2py intent(in) i,c,w,x
      pi=3.1415926535
      gaussian=i/(w*sqrt(2*pi))*exp(-(c-x)**2/(2*w**2))
      end function gaussian
